{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Grid filtering in PET\n\nHow filter work in py eddy tracker. This implementation maybe doesn't respect state art, but ...\n\nWe code a specific filter in order to filter grid with same wavelength at each pixel.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib import pyplot as plt\nfrom numpy import arange\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\n\n\ndef start_axes(title):\n    fig = plt.figure(figsize=(13, 5))\n    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\n    ax.set_aspect(\"equal\")\n    ax.set_title(title)\n    return ax\n\n\ndef update_axes(ax, mappable=None):\n    ax.grid()\n    if mappable:\n        plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "All information will be for regular grid\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = RegularGridDataset(\n    data.get_demo_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"),\n    \"longitude\",\n    \"latitude\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Kernel\nShape of kernel will increase in x, when latitude increase\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 8))\nfor i, latitude in enumerate((15, 35, 55, 75)):\n    k = g.kernel_bessel(latitude, 500, order=3).T\n    ax0 = fig.add_subplot(\n        2,\n        2,\n        i + 1,\n        title=f\"Kernel at {latitude}\u00b0 of latitude\\nfor 1/8\u00b0 grid, shape : {k.shape}\",\n        aspect=\"equal\",\n    )\n    m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap=\"viridis_r\")\nplt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Kernel along latitude\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 8))\nax = fig.add_subplot(\n    111,\n    ylabel=\"Kernel weight\",\n    xlabel=\"Latitude in \u00b0\",\n    title=\"Kernel in latitude, centered at 0\u00b0 of latitude \",\n)\nk = g.kernel_bessel(0, 500, order=3)\nk_lat = k[k.shape[0] // 2 + 1]\nnb = k_lat.shape[0] // 2\nax.plot(\n    arange(-nb * g.xstep, (nb + 0.5) * g.xstep, g.xstep), k_lat, label=\"Bessel kernel\"\n)\n\nax.legend()\nax.grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Kernel applying\nOriginal grid\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will select wavelength of 300 km\n\nLow frequency\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT low frequency\")\ng.copy(\"adt\", \"adt_low_300\")\ng.bessel_low_filter(\"adt_low_300\", 300, order=3)\nm = g.display(ax, \"adt_low_300\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "High frequency\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT high frequency\")\ng.copy(\"adt\", \"adt_high_300\")\ng.bessel_high_filter(\"adt_high_300\", 300, order=3)\nm = g.display(ax, \"adt_high_300\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Clues\nwavelength : 80km\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g.copy(\"adt\", \"adt_high_bessel\")\ng.bessel_high_filter(\"adt_high_bessel\", 80, order=3)\ng.copy(\"adt\", \"adt_low_bessel\")\ng.bessel_low_filter(\"adt_low_bessel\", 80, order=3)\n\narea = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Spectrum\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum\")\nax.set_xlabel(\"km\")\n\nfor label in (\"adt_high_bessel\", \"adt_low_bessel\", \"adt\"):\n    lon_spec, lat_spec = g.spectrum_lonlat(label, area=area)\n    mappable = ax.loglog(*lat_spec, label=f\"lat {label}\")[0]\n    ax.loglog(\n        *lon_spec, label=f\"lon {label}\", color=mappable.get_color(), linestyle=\"--\"\n    )\n\nax.set_xlim(10, 1000)\nax.set_ylim(1e-6, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Spectrum ratio\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum ratio\")\nax.set_xlabel(\"km\")\n\nfor label in (\"adt_high_bessel\", \"adt_low_bessel\"):\n    lon_spec, lat_spec = g.spectrum_lonlat(label, area=area, ref=g, ref_grid_name=\"adt\")\n    mappable = ax.plot(*lat_spec, label=f\"lat {label}\")[0]\n    ax.plot(*lon_spec, label=f\"lon {label}\", color=mappable.get_color(), linestyle=\"--\")\n\nax.set_xlim(10, 1000)\nax.set_ylim(0, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Old filter\nTo do ...\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}